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Abstract 

<""] This paper proposes a novel technique to reduce the computational burden associated with 

the simulation of localised failure. The proposed methodology affords the simulation of damage 

S initiation and propagation whilst concentrating the computational effort where it is most needed, 
i.e. in the localisation zones. To do so, a local/global technique is devised where the global 
(slave) problem (far from the zones undergoing severe damage and cracking) is solved for in 
a reduced space computed by the classical Proper Orthogonal Decomposition, while the local 
^ (master) degrees of freedom (associated with the part of the structure where most of the damage 

f — ■ is taking place) are fully resolved. Both domains are coupled through a local/global technique. 

This method circumvents the difficulties associated with model order reduction for the simulation 
of highly non-linear mechanical failure and offers an alternative or complementary approach to the 
development of multiscale fracture simulators. 
QO Keywords: Adaptive Model Order Reduction (MOR); Local/global Approach, Nonlinear Fracture 

Mechanics; Proper Orthogonal Decomposition (POD); Newton/Krylov Solver; 



1 Introduction 

Simulating damage initiation and subsequent global structural failure is one of the most active topics in 



computational mechanics. Several mathematical models and numerical methods have been developed 
over the years to assess various limit states such as failure due to permanent deformations, cracks 
or decohesion/delamination, e.g. in composite materials. Yet, these models, be they damage based 
or relying on discrete cracks are generally computationally expensive, as they require a fine scale 
description of the structural and material properties. Therefore, today's engineers are not able to use 
these state-of-the-art models for routine design. For important recent advances in the treatment of 
material failure (e.g. discontinuous fracture [T], advanced damage models [2J [3J, damage plasticity 
models [H[5] or their combination [5], etc.) to become useful in practice, it is thus important to devise 
techniques which are able to significantly reduce the computational effort required without sacrificing 
accuracy. 

Historically, reducing the computational time associated with solving nonlinear problems in solid 
mechanics has mainly been addressed by developing homogenisation techniques [3|S]. In this case, 
the material properties associated with a material point in a coarse representation of the structure 
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is obtained by averaging of the fine scale material behaviour over a "representative volume element" 
which suppresses the need to resolve the fine scale explicitly over the (possibly large) structure. These 
methods are widely adopted (see e.g. for a recent description [5]) and are rather effective in decreasing 
the computational costs associated with solving fine scale problems. 

However, the application of homogenisation-based methods to failure simulations is not straight- 
forward, as the assumptions required to prove the separation of scales are not satisfied. This may 
explain why such methods are still under-developed for fracture problems. The ability to coarse-grain 
fracture information would however lead to a further decrease of computational time by concentrating 
the computational effort where it is needed most. In the wake of a propagating crack, for example, 
the full resolution used to resolve fine scale effects close to the crack front is no longer necessary, and 
a coarse description may be sufficient. This issue is being intensively investigated by several teams, 
namely [TO] who developped a coarse graining method called Multiscale Aggregating Discontinuities, 
and [TJJ [T2J H3] who were able to derive an homogenisation scheme for lumping the fine scale deffects 
of Representative Volume Elements (RVEs) into cohesive cracks introduced at the coarse scale. The 
solutions proposed so far are however highly problem-dependent and usually do not inherit the solid 
mathematical foundations from homogenisation theory. 

We propose here to follow an alternative route which relies on very recent advances in model order 
reduction strategies for highly nonlinear structural problems [21 [TSl HH [T71 [18]. We believe that, 
though rarely studied in the context of failure simulation in nonlinear mechanics, this avenue provides 
a possibly interesting alternative (or complementary) approach to homogenisation-based solution pro- 
cedures. This paper is in substance a significant improvement of the work proposed in [18], dedicated 
to the fast prediction of the onset of structural failure by adaptive model order reduction. We aim here 
at preserving the computational efficiency and accuracy obtained in [18] when the size of the damaged 
zone becomes significant compared to the overall structure. 

The definition of model order reduction can vary depending on authors. We speak here of an 
approximation of the unknown field in a coarse space spanned by a basis of global vectors (as opposed 
to locally supported shape functions, in finite elements for example). Model order reduction (MOR) 
finds its roots in the systems engineering community, where the issue of obtaining a fast and reliable 
approximation of transfer functions has been extensively addressed over the years. The interested 
reader can refer to the review proposed in [JJ5]. In the context of solid mechanics, most studies have 
focused on the reduction of dynamics problems. Starting from classical modal basis [20] in a Ritz- 
Galerkin framework, two main classes of model order reduction methods have been proposed. The 
first propose a truncation of the modal basis and an approximation of the high-frequency behaviour by 
static corrections (see for instance [2JJ). The second use the idea of domain partitioning. Slave degrees 
of freedom are eliminated by Guyan reduction (or static condensation) . The famous Craig-Bampton 
method [22] is an evolution of this idea, using truncated modal bases to enhance static condensation. 
These techniques have been extended to the approximation of nonlinear dynamical behaviour (an 
interesting review is proposed in |23j). 

The literature concerning the application of model order reduction to the approximation of the 
(quasi-)static behaviour of nonlinear structures is a lot sparser. Different routes have been investigated 
to develop such reduced models. The first family of methods consists in taking advantage of the a 
priori knowledge of a set of representative solutions to fasten the solution process. These representative 
responses, called snapshots [23], can be for instance a set of solution vectors obtained for given values of 
structural parameters. The Proper Orthogonal Decomposition (POD) [25] [2S] can classically be used 
to extract a relevant global basis from these snapshot vectors, which in turn defines a coarse space 
where the solution to the structural problem with a different set of parameters can be cheaply searched. 
The second family is not exactly based on reduction by projection. Spectral approaches referred to 
as Proper Generalized Decomposition (PGD) are a family of solvers capable of building the solution 
to a multidimensional reference problem as a sum of products of functions of separate variables. The 
dimensions can be space or time, but also structural parameters. Once this solution is obtained offline, 
the online computation simply consists in the evaluation of the separated functions, which are known 
explicitly. These novel methods are very appealing and were recently applied to linear, stochastic 
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and/or time dependant problems [UJ US] (including linear predictions of nonlinear solvers and 
to non- linear partial differential equations [3U1 HE] • Their extension to irreversible constitutive laws 
such as the ones considered in damage mechanics does not seem to be straightforward, and, as far as 
we know, has never been adressed before (it is therefore a very interesting possible extension of the 
present piece of work). 

We focus here on the first family of solvers. To the authors' knowledge, the application of model 
order reduction methods to failure simulations has virtually never been addressed. The main reason 
for this is that the solution fields in the zones of the structure where damage or failure initiate and 
propagate are highly parameter-dependent, and cannot be known in advance. Hence, precomputed 
reduced bases are unable to represent the structural behaviour beyond the onset of strain localisation, 
damage and fracture. In the particular case of linear elastic fracture mechanics, the authors of [311 have 
proposed a mesh morphing technique that allows to keep an autosimilar description of the singularity 
during propagation. In this particular case, the use of a snapshot POD is relevant and allow for 
interesting computational savings. Corrective POD algorithms have also been proposed to alleviate 
this in the case of plasticity or damage models [TSJ [18] . Both these techniques are based on global 
corrections of the initial reduced basis by Krylov algorithms. Though efficient for global nonlinearities, 
it was shown in |18j that in the case of damage assessment, the number of corrections increases when 
nonlinearities localise, i.e.: when failure occurs. The purpose of this paper is to circumvent this 
phenomenon. 

The first stage of the proposed method is to only reduce the set of balance equations which exhibits 
a "smooth" nonlinearity. In the context of the simulation of damage, one part of the degrees of freedom, 
corresponding to highly damaged zones, will be considered as master degrees of freedom, which will 
be fully solved for, while the remaining will be approximated as a linear combination of reduced basis 
vectors, which may be precomputed and corrected on the fly. In this respect, the proposed method 
has strong links with Craig-Bampton methods, or model order reduction by substructure such as the 
one proposed in [32l [33l [29] . A similar approach was also proposed in [34] . The second stage of our 
developments consists in developing an efficient solver for the linearized system of equations obtained 
by a Newton solution algorithm for the reduced problem. The system is condensed on the master 
degrees of freedom, and solved by a Krylov algorithm. The global reduced basis is used to efficiently 
solve the consensed problem: an augmentation technique ensures that the iterations are performed in a 
space which is orthogonal to the one spanned by the reduced basis. As a result, the augmented Krylov 
iterations correspond to a two-level solver, the coarse level being the space spanned by the restriction 
of reduced basis vectors to the fully resolved set of degrees of freedom. 

The paper is organized as follows. In the first section, we recall the basics of projection-based model 
order reduction applied to the solution of nonlinear quasi-static problems. In section 3, we formally 
describe the proposed local/global reduction algorithm. This technique is applied to the reduction 
of damage models in section 4. Finally, we show in the last section that this technique permits to 
drastically reduce the number of corrections of the reduced basis required to obtain a desired level of 
accuracy. 

2 Model order reduction by projection: general concepts 

Projection-based model order reduction can be introduced as a means to obtain, at low computational 
cost, a good approximation to a set of nonlinear equations of the form: 



Where F is a vectorial, possibly nonlinear, function of a set of vectorial state variables X. The 
expression of this function can take various forms, depending on the application. 




(1) 



3 



2.1 Problem statement in nonlinear continuum mechanics 



We consider a structure occupying a continuous domain Q with boundary <9S1. It is subjected to pre- 
scribed displacements C/ D on its boundary dQ u and to prescribed tractions F D on the complementary 
boundary dQf = dfl\d£l u , over time interval [0, T]. Let u be the unknown displacement field, which 
belongs to the space U of kinematically admissible fields: 

U = {ueH\n) \u ldUu =u D ) (2) 

Let U° be the associated vector space. 

U° = {a€H 1 (n) I u [aQu =0} (3) 

Under the assumptions of quasi-static evolution of the structure and small perturbations, the weak 
form of the equilibrium and constitutive laws read, at any time t G [0, T]: 

Vw* G U°, find m€W such that: 



/ g ■. g(u*) dn= [ i u* dn+ I f u .u* dr 



(4) 



'T<t 



where q_ is the Cauchy stress tensor and e(u) is the symmetric part of the displacement gradient. The 
constitutive relation between a and e(u) is nonlinear and described using internal variables (plasticity, 
damage for instance). It is assumed local and rate-independent. 

2.1.1 Space and time discretization 

Let us perform a standard finite element approximation of the space of unknown displacement field hi 
(and a similar approximation of the space of test functions hi ): 

U h (n) = lu(x) | u(x) =Y^Ni(x)u\ (5) 

where n n is the number of nodes, x is the position vector, (^Vj)ie[i,n„] are the finite element shape 
functions associated to each node of the finite element mesh, and (Mi)ie[i,n„] are the nodal values of 
the displacement field. 

The nonlinear solution strategy used to solve the problem over time is a classical time discretization 
scheme for quasi-static and rate-independent problems. This procedure consists in finding a set of 
consecutive solutions at times (tra)ne[o,n t ]- 

The introduction of the finite element approximation and time discretization into equation Q at 
any time t n+ \ of the analysis leads to the following nonlinear vectorial equation: 



Fmt AU, U |tm l+F Ext =0 (6) 

where the vector of increment in the nodal unknowns AU G M" u (n u is the number of scalar nodal 
unknowns) is defined by AU = Uu — U.u , F Int G K n " and F Ext G M n " are respectively the internal 
forces resulting from the discretization of the internal virtual work (left-hand side of the first equation 
of system Q) and the external forces resulting from the discretization of the external virtual work 
(right-hand side of the first equation of system Q). In the following, the dependency of the internal 
forces with respect to the history of the unknown fields will not be written explicitly, unless necessary. 

The set of successive solution vectors to problem ^ at times (£n)ne[o,n t ] i obtained by a tangent 
Newton algorithm, will be considered as the reference solution. The algorithms used in the following 
should provide a solution that is closer to the reference solution, but for a cheaper computational cost. 
Therefore, the relevance of (i) the initial space and time discretizations and (ii) the chosen Newton 
solution algorithm will not be discussed in this paper. 
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2.2 Reduced solution of the balance equations by projection 
2.2.1 Principle 

The solution vector is searched for in a space of small dimension (several orders smaller than the 
number of finite element degrees of freedom). Let us call C the matrix whose columns form a basis of 
this space: 

Q = ( c 1 C 2 ... C n " ) (7) 

where n c is the dimension of the reduced space, and (C fc )/ cS ji w j € (M n ")™ c are the basis vectors. 
Applied to the reduction of problem ([6]), the increment in the solution field is approximated by: 

AU=£a (8) 

where we introduced the reduced state variables a 6 R™ c . Problem ^ is now overconstrained. One 
usually tackle this problem by prescribing a Galerkin orthogonality condition: the residual of equation 
([6| R = F Int (U| tn + C a) + F Ext is constrained to be orthogonal to any test vector SXJ * = CSa* 

belonging to the space spanned by the reduced basis vectors (C )fc e [i i7lc j '■ 



The reduced form of problem Q is: 



C T R = (9) 



(F Ext +F Int (U| tn +ga)) =0 (10) 



2.2.2 One particular choice of projection space: the proper orthogonal decomposition 
(POD) 

Various choices of projection spaces are proposed in the literature. One of the most successful, as far as 
nonlinear simulations are concerned, is the snapshot-POD [24. This particular technique requires the 
knowledge of a representative family of solutions to the global problem. This set of vectors (S fc )fcgji jns ] 
is called snapshot. The aim is to find an orthonormal basis {C k )kefi,n c j > of dimension n c smaller than 
n s such that the distance between spaces Im(C) and Im(S) is minimum. This distance can be quantified 
by the following metric: 



3 = 1 



(11) 



which must be minimized under the constraint of orthonormality of (C )fee[i,n e |- Problem (111 can 
be reformulated as an unconstrained formulation: 

L(C\ . . . , C n % An, ■ , A„ c „J = J(C\ . . . , C n ") + £ £ (c iT C j - Sij) (12) 

i=l j=i 

where (Ay)j i;J -g[i n j are Lagrange multipliers and Sij the Kroncckcr symbol. The n c (n c + l)/2 optimal- 
ity conditions for maximizing the Lagrangian with respect to the lagrange multiplier naturally ensure 
the orthogonality of the basis: 

FIT t 

^ = -> C"'C- i tJ (13) 

Whereas the n c optimality conditions for minimizing the Lagrangian with respect to the basis vectors 
reads: 

|^ = -> (c lT sA = X vl C l and A„ =0 for i^j (14) 

d ± J= i 
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Defining A, = X u and § = ( S 1 S 2 ... S" s ) € R n " x E n », equations ((T4|) can be written as the 
following n u x n u eigenvalue problem: 

gg T C l = X, C l (15) 

This problem is equivalent to computing the Singular Value Decomposition (SVD) of S, where the 
singular values Si are such that sf = Xi > 0. If n s <C n u (which is the case in general when applying 



the snapshot-POD), replacing problem (15 1 by the eigenvalue problem on the n s x n s operator S S is 
computationally cheaper, and reads: 

C^A^IY 1 with I^V^AiV* (16) 

The error associated with the truncation of the SVD at rank n c is quantifyed by the distance J(C 1 , C™ c ) 
which can be shown to be equal to the sum of the truncated eigenvalues: 

n s 

J(C\...,C n °)= £ Xi (17) 

i—n c + l 

This result provides a precious tool for choosing the most relevant basis. Indeed, only the basis vectors 
Cj associated with the largest Xi (in practice, those such that Xi/X ma x > s, e being a small parameter 
depending on the required accuracy) are selected as reduced basis vectors. This ensures that for an 
given number n c of basis vectors, the approximation of the snapshot space is maximised in the sense 
of p|. 



2.3 Nonlinear solution algorithm 



At any time step, problem ( 10 1 can be solved by classical Newton algorithms (successive linearizations). 
Iteration i + 1 of the tangent Newton method consists in solving the following linear predictor: 



(18) 



where the increment in the reduced state variables is Sa l+1 = a. i+1 ~a\ and the residual and tangent 
operator obtained from the knowledge of the solution at iteration i read (see [18] for more details) : 



F Int (U |t „ +£«*) + 
i„t(U| t „ +Qa) 



-Int 

<9F 



-Ext 



(19) 



=T,R 



da 



3 Local/global model order reduction strategy 

Projection-based model order reduction techniques introduce a global approximation of the displace- 
ment. As such, even adaptive versions (35l [T5] of these methods are not well suited to the analysis 
of localised nonlinearities. We propose in the following to use reduction techniques by projection for 
the "weakly nonlinear" equations of reference problem ([6]), namely equations for which reduced order 
modelling (ROM) is relevant, while the remaining equations will be solved directly. 



3.1 Displacement approximation 

The principle of the proposed local/global strategy is to split the unknown solution vector into two 
parts, only one of them being approximated as a linear combination of reduced basis vectors. We 
shall use superscript ^ for the approximated "slave" part of the solution vector, while superscript ^' 
(which stands for "fully resolved") will be used to denote its complementary "master" part. 



G 



More precisely, let AU be the unknown increment vector at time t n +i, the numbering being 
reorganised as follows: 

m= [m^) = {^)— (20) 

O) — (/) 

In the above notations, AU G M™ r and AU G M. n f , with n r + rif = n u and rif much smaller than 
n r . § (r) G {0, l}"- x {0, 1}"" and ]| (/) G {0, 1}"' x {0, l} n * are two boolean extractors (rectangular 
matrices with one 1 per line, the other coefficients being null). 

AU is approximated as a linear combination of global vectors (C fc )fce[i,n c ] G (]R n ")™ c : 

AU (r) = (0 r) Q) a (21) 

where a G M™ c is the vector of reduced degrees of freedom. In the initial numbering, the unknown 
vector can be expressed in the form: 



AU = AU (r) + AU (/) 

where 



AU (r) = E (r)T AU (r) (22) 



v AU (/) = i (/)T AU (/) 

Introducing projector = E*- r ^ E^ r ^ , the approximation of the displacement increment finally reads: 

AU = P (r) C a + g( / ) T AU (/) (23) 
For the sake of legibility, we define the new state vector X at time t n +\- 

s-(~b) < 2 ") 

Using the above definition, the approximation of the displacement increment is now: 

AU = AX where A = (P (r) C E (/)T ) (25) 



3.2 Locally reduced set of balance equations 

At that point, we need to define how the overconstrained balance equations ^ will be solved for. Let 
us restrict our numerical investigations to the Galerkin procedure. The balance equation are therefore 
required to be orthogonal to any test vector: 

<*IT = A ( —(/)* j (26) 

* — (/)* 

where 5a G K" c and (5U G M n/ are arbitrary test state variables. The orthogonality requirement 
applied to the initial nonlinear set of equations ([6]) yields the following nonlinear system: 

R R (X) = A T (P Int (aU(X) + U |t J + F Ext ) = (27) 
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3.3 Newton solution scheme 



At each time step of the time discretization scheme, the reduced problem (27) is solved by a Newton 
algorithm. The (i + l) th Newton iteration consists in finding SjC +i 
n c + m linearized equation: 



0Er(X) 



ax 



R 



R 



satisfying the following set of 

(28) 



X=X ! 



where <DC +i = X ,+1 — X' an d Sir is the residual of problem ([27]) computed at iteration i: 



The expression of the tangent operator K 



Sr — Rr(X 4 ) 

5R R (X) 



(29) 



T,R 



ax 



is obtained by differentiation of 



X' 



equation (27) with respect to the set of reduced state variables, in any arbitrary direction <5X such 
that: 

(«) 



This differentiation reads: 



9R r (x) 



dX 



8X 



Using the tangent operator of the full system of equations ([6]), and assuming that 
write: 

, T 9F Int (AU) <9AU 



dF 



Ext 



<9AU 



i;,R^ = A 



dAU 



AU=AX ! 



ax 



<5X 



(31) 
0, on can 
(32) 



AU=AX' 



Using the notation K 



dF 



Int. 



the previous equation yields the expression of the reduced 



AU (X') 



tangent operator as a function of the tangent of the full system of equations: 



(33) 



Introducing the above expression of the tangent into equation (28), and expanding the expression of 
operator A, one finally obtains the following linearized system: 



E (/) K^P (r) C E (/) K^E (/)T 



i+i 




SX 



Fi„t(AU(X l )+U |t 



(34) 



-Ext 



K" e 

=T,R 



Notice that problem (34) is considerably reduced compared to a direct linearisation of Indeed, 
= x W lf where n c <C n u and, applied to the analysis of localised nonlinearities, we can 
reasonably expect that n/ is at least one order of magnitude smaller than n u , which means that only 
a few of the balance equations are strongly nonlinear. This assumption will be validated in section [4] 
This framework is very general. It does not rely on domain decomposition but on a splitting of 
the balance equations. Yet, this technique will be particularized in section [4j and we will show that 
using the ideas inspired from domain decomposition methods 36, 37 j and local/global methods 38, 39J 
provides efficient splitting in the case of damage mechanics. 
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3.4 Two-level solution of the linearized systems 



The successive solution to ( 28 ) could be obtained by a direct solver. This would not affect the results 
given in the following. However, we believe that, for this particular application, an iterative solver has 
several advantages: 

• Iterative solvers, such as a conjugate gradient (see [10] for instance), do not require any fac- 
torisation. Factorisations usually deteriorate the sparsity of algebraic systems, which results in 
an increase in the memory requirements. Newton-Krylov methods do not, in theory [41] . even 
require the assembly of the stiffness operator. 



Direct solvers are less versatile than iterative ones, and are usually a limiting point to extend 
numerical solution schemes to parallel computing. 



Solving exactly the successive linearized systems (281 is not required. This idea has led to 
intensive work on Inexact Newton Methods [42 , but this issue will not be addressed in this 
paper. The interested reader can refer to [37J HH] for studies on similar topics in the scope of 
domain decomposition and model order reduction. 

More importantly, a basic direct solver such as Cholesky factorizations will not make use of any 
mechanically relevant information obtained during the previous iterations of the nonlinear solver. 
In the chosen iterative algorithm, the part of the displacement field that is "fully resolved" will 
be solved for efficiently using the reduced basis as a preconditioner, following the ideas proposed 
for instance in [43 ] [4T ] [37] . 



We use here a Krylov algorithm to efficiently find a solution to linearized system of equation ( 28 ) 
Three stages of preconditioning are performed: a condensation, a projection and a crude precondition 
ing of the resulting system. 

3.4.1 Condensation 



Let us condense the linearized problem (34) on the master (i.e.: "fully resolved") degrees of fredom: 



i^5U (/) =R^ (35) 

where the superscript corresponding to the current and previous Newton iterations have been dropped 

is the primal Schur complement and Rg ' 



for the sake of legibility. is the primal Schur complement and R^ the condensed residual. They 
are defined by 

Where we have used the usual block notations: 



/) 

- T ' R (36) 



=T,R =1 
K (M K (, 

=T,R =T,R 



E R = ( J (38) 



The state variables corresponding to the reduced part of the displacement vector increment have 
been eliminated from the system of equations and can be retrieved as follows: 

6* = ~ (Kg) 1 + K£2 (39) 
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(^"T r) ^ s an °P era ^ or °f ver Y small size, which can be explicitly computed prior to applying 
the iterative Krylov solver. However, in the case where a very large number of reduced basis vectors is 
used to approximate the unknown field in the smooth region, computing the Schur-complement/vector 
products corresponding to the construction of the successive Krylov basis vectors will become compu- 
tationally more efficient. 

3.4.2 Iterative solution to the condensed problem 

The classical projected (or augmented) conjugate gradient [44j [36] is applied to the approximate 



solution of (35) (the linearized operator is assumed symmetric, positive and definite). The chosen 
augmentation space is the space spanned by the reduced basis vectors Im(C^) with the restriction 

The starting point of augmented Krylov solvers is to decompose the unknown solution increment 
into two supplementary spaces: 



where 



SVc ) € Im (40) 
5U« e Im(£^) V> = Ker (fi^ T §^) 

lo(/j designing the Sj^' —orthogonality, which is ensured by introducing the classical projector: 

E = Id - CW(C^ T S^CW)- 1 cW) r |(/) (41) 

and writing that 8U$ = PgUjjft. 

This separation of the search space into two subspaces Im fc^J and Im(P) in direct sum leads 
to the following uncoupled equations: 



(p T i? ) p)«4 /) =g T s^ ) 



(42) 



The first line of equation (42 1 is a coarse initialization of the projected conjugate gradient. The second 

line is the linear prediction of the full problem projected on Im(C^) =p . This system is symmetric 
and can be solved by a preconditioned conjugate gradient. Hence one solves iteratively: 

M _1 (if'ip 3 P) §H { k = g T Ec /) (43) 

where M is the left preconditioner (symmetric, definite and positive). In our test cases, M is a 

diagonal matrix whose entries are the elements of the diagonal of S j/' . Algorithm jlj gives the classical 
implementation of an augmented preconditioned conjugate gradient. 

3.4.3 Interpretation and comments on the expected additional costs 

The three preconditioning techniques that we have detailed in the previous sections define a cheap 
structural preconditioner. A schematic representation of this algorithm is given in figure ([!]). 

From an algebraic point of view, it can easily be proved that the condensation has the same effect 



on the conjugate gradient iterations as a left block preconditioner for system (34) 



(s£2) 











o / 



Mp , 1 = I V=t,r; = ) (44) 
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A priori reduced basis 



New problem 



Coarse initialisation in the reduced space 



iUi 6 Im(E (/) C 



Condensation of the linearized reduced 
problem on the master degrees of freedom 



AAA 



AAA 
AAA 



AAA, 



Computation of the displacement 
of the master degrees of freedom 



-SStP 




1 



Computation of the displacement 
of the slave degrees of freedom 



<5U = E (/) <5U (/) + P (r) C<5a 




AAA 



Orthogonalization and 
Minimisation 



K,j 



K,j-i T oy w, 



Computation of the residual BcGi of the 
projected problem ("non-corrected" search 
direction) 



AAA 



Coarse scale correction of the 
search direction by orthogonality 
with respect to the reduced space 



Z, = PM l R r 





] 



Figure 1: Graphical representation of the two- level solution algorithm. 
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Algorithm 1: Augmented preconditioned conjugate gradient for the solution to the condensed 
linearized reduced problem 



Compute < 



tug 



aV) - v(ff) _ KV r > ( K<< rr A 1 K^ r - 

=P =T,R =T,R \==T,R / =T, 
(fr) fr/Hr 1 R (r) 
T,R \=T,R / — R 



{ Q(f) T gU)QU)y 
1 (/)(C(/) T sj/ ) C (/) ^ 



-i c (/)^ R (/) 



=EiCG,0 — fic 



C/) -S(/>5U^=P T R (/) 



z = PM Rcg,o 



attC/) — n. 



, W = Z and 
for j = l,...,n do 

- (Scgj-lWj-iJ/^Wm.Wh 
aj-iWj-i 



<5U (/) - <5U (/) 



R 



GGj — ±±CG 
Zj = EM^RcGj 



end 



From a mechanical point of view, condensing the problem on the master degrees of freedom ensures that 
any solution provided by the iterative solver exactly solves the reduced part of the balance equations. 

The augmentation of the conjugate gradient can also be seen as a left preconditioner for the 
condensed system of equations (second line of system (|42|), which is associated with an initialization 



of the iterative alogrithm. The initialization S\J ^ (first line of system (42)) is the best solution of 
the condensed problem in the restriction of the pre-computed reduced basis to the master degrees 
of freedom (in the sense of a (Sj/^) _1 -norm). Taking into account the previous comments on the 
interpretation of the condensation step, this initialisation is exactly the solution to the linearized 
system of equation ( 18 ) on the master degrees of freedom. In other words, the initialization of the 
iterative algorithm provides the solution of the reduced linearized problem which would have been 
obtained without local/global enhancement. No significant additional costs compared to the basic 
projection-based reduced order modelling are involved so far. In the case where the solution which is 
looked for is correctly approximated in the reduced space (e.g.: "smooth" nonlinearity, or unloading 
of the structure), no conjugate gradient iteration is required. 

If conjugate gradient iterations are required, the augmentation acts as a global correction of the 
successive search directions (see figure dTJ). It ensures that any solution of the form (40 1 still satis- 
fies the reduced balance equations without local/global enhancement. Therefore, the complementary 
solution SXJ yJ can be interpreted as a local correction to the solution provided by the pre-computed 
reduced model. We will see in the examples that, as the conjugate gradient does not need to search for 
the part of the solution that is directly captured in the pre-computed reduced space, the convergence 
rate is very high. 

The diagonal preconditioner M classically accounts for local heterogeneities in the structure. A 
much more sophisticated preconditioner can of course be used (e.g.:incomplete Cholesky factorization). 
Yet, if the reduced basis is correctly pre-computed and/or updated, the computational effort required 
to obtain a good solution of the projected problem should be relatively low. Hence, the potentially 
expensive construction of an advanced preconditioner is not a priori justified. The example section 
will show that the number of conjugate gradient iterations required to compute a correction to the 
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reduced model is indeed kept very low, even while only using this crude preconditioner. 



4 Application to the adaptive reduction of damage problems 
4.1 Damageable lattice model 

We focus on a very simple damageable lattice structure, made of bars under traction or compression 
(figure |3j)). Elementary bar indiced b occupies domain Of, such that fl — U&e[i n 6 ] ^ fc > wnere n b 
is the number of bars. The direction of bar number b is denoted by n kl = PkPi/\\PkPi\\ where P k 
and Pi (k < I and (k, I) G Jl,np] 2 ) are the two extremities of the bar as defined in figure The 
displacement field in is supposed constant in a cross-section. Therefore, the displacement of any point 
M G fib can be expressed in the following form: 

M|Mefi 6 = (Mffl-tkki) Thti = u (x) m ( 45 ) 

where each cross section is parametrised by a scalar x G [0, ||PfcP;||], defined by x = P}.M.n kl . 

We use a constitutive law based on classical damage mechanics to describe the behaviour of the 
lattice structure. The lineic strain energy of bar b reads: 

e dls = ±E(l-d)S b ef £ (46) 

where eu = u^ix) is the strain measured in the direction of the bar, {Sb)beli,n b } , is the section of the 
bar, and d is a damage variable which ranges from (undamaged material) to 1 (completely damaged 
material). The local state equations, derived from the strain energy read: 



de 



N ls = -^ = E(l-d)S b e l£ 



where we have introduced the axial force N = Sb fajy-<I-?l/ ! z)j and the thermodynamic force Y. A 
non-local evolution law is defined to link the damage variables to the thermodynamic forces: 

d\,f = max I min f 1, max (a(Y\f } } (48) 

We suppose that the lattice is not subjected to any volume force, that the material of each elemen- 
tary bar is isotropic and homogeneous, and that the section of each bar is constant. As a consequence, 
the resulting displacement field evolves linearly along the direction of each bar. 

In terms of finite element discretization, each bar will be considered as a linear element 



Because of the assumptions made previously on the loading, geometry and material properties, the 
finite element solution yields the "exact" solution to the damageable lattice problem. 



4.2 Lattice structure and loading 

We solve the problem defined in figures (J3j) . The sections and Young's modulus are unitary. The 
material parameters are /3 = 0.5 and a = y2. Any bar parallel to axis (0,x), (0,y) or {0,z) has a 
unitary length. 
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Figure 2: Definition of the lattice problem. 




Damage 



0.5 1 

Figure 3: Reference solution obtained at the eighth load increment of the full-scale simulation (max- 
imum value of prescribed forces before unstable damage propagation). Darker bars have a higher 
damage state. 
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The Dirichlct boundary conditions applied to the structure read, for any node (P,) (with i S [1, n p J) 
of the lattice such that (P,) € (0,x,y): 

u {Pi .z = (50) 

As represented in figure homogenous Neuman boundary conditions, in the — z direction, are applied 
to any point Pj that belongs to plane V defined by: 

V = {M = (x, y,z)\x€ [7, 9] and y € [8, 10] and z = 11} (51) 



Damage 




Figure 4: Final damage obained at the end of the full-scale simulation, performed in 30 load increments. 
The bars that are completely damaged are not represented. 

Increasing the value of the applied forces leads to the initiation and propagation of a damaged zone 
as shown in figure Q. This phenomenon is unstable from a structural point of view, and the damage 
state represented in this figure is past the limit point at which the force applied to the structure 
reaches its maximum value. An arc-length algorithm with local control is used to solve the damage 
propagation beyond this limit point. The basic idea of this algorithm is to look for the amplitude of 
the force under the constraint that the maximum increment in the damage variables over a time step 
is fixed (see [351 H] or [IE] for additional details) . The reference solution is obtained by solving the full 
problem ([6| in 30 of these damage increments of the Newton/arc- length algorithm. The final state is 
represented in figure Q. 

Our purpose is to check the adequacy of the proposed local/global reduction technique when trying 
to obtain at cheap costs the load/deflection curve (pre and post failure phases) of the structure under 
the loading conditions described previously. The results provided in the following sections give the 
value of the norm of the external forces vector (also called loading factor) obtained by the arc-length 
method. The displacement used to plot the load/deflection curves is the average vertical displacement 
of the points at which forces are applied. 
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Figure 5: Splitting obtained at two different stages of the simulation performed by the local/global 
reduction method (left: load increment 4, right: load increment 30). Dark bars are connected to at 
least one fully resolved node. As time increases, the degrees of freedom for which no reduction is 
performed localises in the region of failure. 

4.3 Particularisation of the local/global approach 
4.3.1 Definition of the Local/Global splitting 

We need at that point to particularize the general method developed in section [3j The splitting 

(r) (/) 

between reduced degrees of freedom AU and the complementary ones AU is made by considering 
that the part of the structure undergoing strong damage variations will not be correctly solved for 
when projected on a pre-computed global basis (see |38[ I37j for similar assumptions in a domain 
decomposition-based local/global approaches). This statement is purely based on mechanical intuition 
and will be validated in the following examples. 

The following procedure is adopted. The splitting will be fixed over a time step, to avoid the 
stagnation or divergence of the Newton process. At the end of a time step, the element undergoing 
the maximum damage increment is spotted, and a sphere of radius p s , centred at the isobarycenter of 
the element is created. Every degree of freedom belonging to a node located inside this sphere is set 
as a fully resolved degree of freedom. 

This procedure is repeated on the remaining elements (those which do not have all their degrees of 
freedom in the previous list). 

The algorithm is stopped if one of the two following statements is satisfied. 

• The maximum damage increment in the remaining elements is lower than fcoam times the global 
maximum. In our examples, /cDam is set to 0.5. This ensures that the fully resolved degrees 
of freedom are actually connected to elements undergoing a significant increase in their damage 
state compared to the remaining of the structure. 

• The number of fully resolved degrees of freedom exceeds fcLocGio times the total number of 
unknowns. /clocGIo is set to 0.1 in our examples. This statement will be satisfied if the former 
one is not, which means that the structure undergoes a virtually homogeneous increase in its 
damage state. In that case, the local/global procedure might be inefficient. 

To illustrate this procedure, two different reduced domains are represented in figure (|5}. The dark 
elements are connected to at least one fully resolved node. The first one is obtained at the end of the 
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third time step and is used to solve the fourth time step of the analysis, while the second one is used to 
solve the thirtieth and last time step. The first figure is obtained during the initiation phase. Notice 
that the maximum damage increments are found in the zone where the the load is applied, and in the 
"pillar" that is closer to this zone, which undergoes a significant bending loading. In the later case, 
corresponding to unstable propagation phase, the damage increment localises in the zone where the 
load is applied. 

This empirical choice for the local/global splitting is only justified by its efficiency. However we only 
considered one type of model for brittle fracture. We expect the optimal splitting to depend on the 
chosen model and on the stage of the structural failure. A more suitable criterion should be investigated 
in future work, based on algebra to circumvent this dependency. For instance, one could consider the 
decrease rate of the singular values of the correlation matrix restricted to an assumed smoothed domain 
and loosely maximise it, with a greedy-like algorithm, by successive small modifications of this domain. 



4.3.2 Modification of the reduced basis 

At the end of a time step, the reduced basis is enriched in order to take into account the information 
provided by the local iterative solver. The procedure used here is simply to add the successive solutions 
to the initial snapshot, and sort them by a singular value computation. This step can be optimized, 
following our previous work [18] . The issue here is that the reduced basis is only locally enhanced, 
which might results in ill-posed reduced problems if not correctly addressed. This question will not be 
further discussed in this paper, which focuses on the benefits provided by the local/global approach. 



4.4 A posteriori reduced basis 

4.4.1 Reference solution to construct the a posteriori reduced basis 

The first test consists in using as a snapshot the consecutive reference solutions obtained by a Newton/arc- 
length solution algorithm performed on the initial system of equations, without reduction. In other 
terms, we try to reproduce the snapshot by reduced order modelling. It terms of realistic applications, 
this example is not of particular interest as one would need to know the solution in advance in order 
to obtain the followinf results. However, studying this academic case is interesting to understand the 
behaviour of the proposed strategy. 

The snapshot is composed of n s — 30 vectors. This information is compressed by using a singular 
value decomposition as detailed in section 2.2.2 We plot, in figure ([6|, the normalised error in the 
snapshot representation 



s j - 2 (c iT s j ') c' 



2\ 2 



^SVD 




(52) 



as a function of the order of truncation of the SVD, where the eigenvalues (^i)ieli,n s l °f the correlation 
matrix are now sorted in decreasing order. The decrease rate of this criterion is relatively low, which 
suggests that the information from the snapshot is poorly suited to compression. For reference, in 
|15| . the author estimates that the snapshot space is correctly approximated by the space spanned 
by the first n c singular vectors if ^svd < 10 -8 . In our case, this would lead to the selection of the 
30 singular vectors to define a relevant reduced model. The singular vector basis is here truncated 
at order n c — 3 in order to yield a significant (thus observable) error in the solution when using the 
basic POD method. This error can then be compared to the one obtained when using the proposed 
local/global scheme. 
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Figure 6: Measure of the error in the snapshot representation as a function of the error of truncation of 
the singular value decomposition. In solid lines dark line, the set of 30 reference solutions are defined 
as the snapshot. In grey lines, a second snapshot operator is defined by concatenating the 30 reference 
solutions and the 30 solutions of the nearby problem. Both analyses are repeated when excluding the 
process zone (defined a priori by the set of points {M = (x, y,z)\x € [5, 13] and y £ [6, 13]})) from 
the snapshot definition (i.e.: the snapshot now is the restriction of the solutions to the part of the 
domain which does not undergo significant damage). The resulting error curves are plotted in dashed 
lines. 



In figure ([7]), the reference load/deflection curved is compared to the one obtained when using 
the basic POD on the first hand, and when using the local/global reduction technique on the other 
hand. The peak load obtained with the basic POD is slightly underestimated, by less than 2%, but the 
overall behaviour of the structure is correctly captured. The local/global algorithm permits to further 
increase the accuracy of the approximated solution. The error in the maximum load factor is now less 
than 0.5 %. 



We now analyse the effect of preconditioning the local iterations of the conjugate gradient by the 



projection technique detailed in section 3.4.2 The convergence curves that are plotted in figure M 



left) are obtained at the third Newton iteration of the fiftieth time step of the analysis (i.e.: during 
the damage propagation phase). The norm of the condensed system (35) of equation, normalized by 
the norm of the condensed right-hand term is observed as a function of the number of iterations of 
the conjugate gradient solver. One can see that the number of iteration to convergence is significantly 
decreased when using the augmented Krylov solver. 

An other interesting observation is that the convergence curve decreases in a monotonic manner, 
which means that the error in the global residual keeps decreasing as a function of the cost of the local 
conjugate gradient. This is not, of course, a proof, and this conclusion is probably problem dependent, 
but the same behaviour as been observed in all our numerical experiments. A probable reason for 
this is that the augmentation deflates the spectrum of the condensed operator, removing its highest 
eigenvalues and forcing the Krylov algorithm to its super-convergence phase (one can refer to [46] for a 
proper analysis of the effect of eigenvectors-based augmentation on the convergence of Krylov solvers) . 
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Figure 7: Load/deflection curves obtained one using successively direct numerical simulation, simple 
proper orthogonal decomposition and local/global reduction technique. The snapshot used here is the 
solution of the direct numerical simulation. 



4.4.2 Nearby solution to construct the a posteriori reduced basis 

In the second test case, we use as a snapshot the consecutive solutions obtained by solving a nearby 
problem, which differs from the reference one only by the position of the applied forces: homogenous 
Neumann boundary conditions are applied to any point Pi (with i G [1, npj) such that Pi S V defined 
by: 

V' = {M = (x, y,z)\x£ [9, 11] and y e [9, 11] and z = 11} (53) 

This snapshot is for instance a particular output of a parametric study meant to obtain the response 
of the system subjected to various load cases. One would expect to be able to reuse the information 
generated during this simulation to fasten the solution process of the next one. The difference between 
the snapshot solution and the reference solution that is looked for is illustrated in figure It 
should be noted that far away from the part of the structure where Neumann boundary conditions are 
applied, the solutions are qualitatively similar. Yet, they are very different within this process zone, 



which qualitatively justifies the local/global splitting proposed in section (4.3.1). The same effect 
appears in figure ([6j which shows that the combination of the information from the two load cases can 
be significantly compressed when removing the process zone from the SVD analysis. 



Figure (10) compares the load/deflection curves obtained in the reference case, and when using 
successively the basic POD strategy and the local/global reduction technique. Obviously, the POD 
strategy looses its relevance very quickly, as the damage localises in the part of the structure where 
the damage would localise during the snapshot simulation. The resulting model is too stiff, and the 
maximum strength of the structure is overestimated. When using the local/global approach, with the 
parameters described previsouly for the fine analysis of the damaged areas, a much more accurate 
load/deflection curve is obtained. Yet the peak load is still overestimated by 5%. This remaining 
inaccuracy is due to the fact that the infomation from the snapshot is not correct in the reduced zone. 
Global corrections of the reduced basis are still necessary. This issue is discussed in the next section. 

In figure ([8j right) we analyse the effect of preconditioning the local iterations of the conjugate 
gradient by the projection technique. As in the previous test case, the convergence curves are observed 
at the third Newton Newton iteration of the fiftieth time step of the analysis. Again, the number of 
iteration to convergence is significantly decreased when using the augmentation. We would expect this 
precontionning technique to be less efficient in this particular case. Yet it is not the case, which is 
probably due to the fact that the global update described in section [4. 3. 2| is mechanically relevant. 
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Figure 8: Convergence curves of the augmented conjugate gradient used to solve the local problem, at 
the third Newton iteration of load increment 15. On the left-hand side, the reference solution is used 
as a snapshot. On the right-hand side, the snapshot used to define the reduced model is the solution 
of the nearby problem. 



5 Adaptive reduced basis 



We now try to reduce the inaccuracy observed in figure ( 10 ) by using global corrections of the reduced 
model. The adaptation algorithm used here was introduced in [T8]. We just recall the basics of this 
particular strategy, and explain how it can be used to enhance the proposed local/global scheme. 

5.1 "On-the-fly" global corrections of the reduced basis 

The principle of these global corrections, is to update the reduced basis used to define reduced problem 



(10). It is, for now, not linked to the global/local technique described previously in this paper. The 



updates are done "on-the-fly" , at any iteration i + 1 of the Newton algorithm used to solve ( 10 ) at 
a given time step t n+ \ of the analysis. If the reduced problem is sufficiently converged, and if the 
residual of the norm of the initial problem evaluated at iteration i is estimated too high, then the 
following linearized system is considered: 

Kt<SI = ~R (54) 

K is an approximation of the tangent K^, while R l is the residual of the initial system of equations 

^ evaluated in Uu + C l a 1 . Operator C* designates the reduced basis considered at iteration i of 
the Newton process. 

The solution to this problem is searched for in two supplementary spaces by a projected Krylov 
algorithm: 

SXJ = £U C + <5U K 

f SJoGlmCC*) (55) 
where < ■ ±-i x-i 

\ <5U K G Im(£ ! ) ^ = Kcr(£ I J^) 

This problem is solved coarsely by a projected Krylov solver (10 _1 is the typical value of the stopping 
criterion), othogonally to the current reduced basis. The resulting solution, which is — orthogonal 
to the CT, is added to the reduced basis: 



ti+l _ I c-ii 



(56) 
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Figure 9: Two set of reference solutions obtained by direct numerical simulations and used as snapshot 
vectors in our examples. Only the solutions to the eighth load increment are represented in both cases. 
On the left-hand side, the snapshot is the reference solution which is looked for. On the right-hand 
side, the snapshot is the solution to a slightly different problem ("nearby problem"), where the position 
of the applied load has been modified. 
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Figure 10: Load/deflection curves obtained one using successively direct numerical simulation, simple 
proper orthogonal decomposition and local/global reduction technique. The snapshot used here is the 
solution of the direct numerical simulation to the nearby problem. 
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Newton prediction i + l can now be solved after updating residual R^ and reduced stiffness K^, with 
respect to this new basis vector. 

Various features are added to this algorithm to efficiently sort the successively added basis vectors, 
but the main idea described above is sufficient to understand the following results. 

5.2 Coupling with the local/global approach 



Once the global norm of the residual of locally reduced problem (27 1 is sufficiently converged, and 
if the global residual is juged too high, a correction vector to the global reduced basis is calculated, 
following the algorithm given in the previous subsection. The residual and tangent operator of problem 



(27) are updated and condensed tangent problem (35) is solved with the projected conjugate gradient 



algorithm. This algorithm converges in a very few iterations as the new basis vector, added to the 
augmentation space, is very close to the solution which is looked for. The slight mismatch is only due 
to: 

• the approximation of the tangent used to provide the global correction vector. 



the coarse solution of ( 54 ) 



• the approximation of the displacement in the reduced zone of the local/global approach. 
5.3 Results 



The test case described in section 4.4.2 is solved using successively the local/global technique on its 
own, the global correction algorithm on its own, and the local/global reduction technique coupled 
with global corrections. The target norm of the relative global residual norm is 10 _1 . Corrections are 



allowed only if the relative norm of the residual of locally reduced problem (27) is below 10 . ft is 



clear from the load-displacement curve in Figure 11 that the global corrections significantly improve 



corrections is given in figure ( 13 1 



the results of the local/global algorithm. (11). An illustration of the POD basis and its successive 



More importantly, the number of relatively expensive global corrections required to achieve the 
target in terms of global residual norm is very low when using the local/global approach, as opposed 
to the one observed when using only global corrections. This effect is shown in figure ( |12[ ). On can also 
notice that using only the global correction methods lead to an increase in the computational costs in 
the damage propagation phase of the analysis, which has been reported in [18 . With the local/global 
reduction approach, this effect disappears. The same trends are observed in [37], where a local/global 
algorithm was developed in a domain decomposition framework to avoid unnecessary computations far 
away from the process zones. This suggests that the result obtained here may be extended, and that 
the number of global corrections required to achieve a given level of accuracy does not depend on the 
local nonlinearity (i.e: independent on the damage state). 



6 Conclusions and discussion 

The work described in this paper focuses on decreasing the computational effort required in solving 
large scale damage and fracture problems, where small scale phenomena must be taken into account to 
accurately represent the global structural behaviour. Although most investigations on this issue rely 
on the extension of homogenization frameworks, we followed an alternative route, choosing projection- 
based model order reduction as a starting point. The proposed local/global reduction technique was 
demonstrated on a damaging 3D frame structure, but we expect that the algorithm will carry forward 
to the more general case of the initiation and propagation of cracks in structures. 
The main conclusions of the paper are: 
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Figure 11: Load/deflection curves obtained one using successively direct numerical simulation, lo- 
cal/global reduction technique and local/global reduction technique with "on-the-fly" global correc- 
tions. 
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Figure 12: Number of global corrections performed when using the adaptive model order reduction on 
the one hand, and when using adapative local/global model order reduction on the other hand. 
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Figure 13: A posteriori POD basis (3 vectors) and successive global corrections (4 vectors computed 
respectively at time steps 1, 3, 10 and 18). Two different views of each basis vector are represented 
for better qualitative understanding. 
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• The proposed local/global model order reduction approach can significantly improve the rele- 
vancy of using a global reduced basis to approximate the displacement field in the case of localised 
failure; 

• The number of global corrections of the reduced model required to obtain a given level of accuracy 
is drastically reduced when the balance equations which exhibit the highest level of nonlinearity 
are excluded from the reduction. 

In most of the applications for which this methodology is intended, the additional computational 
costs are kept minimal for two reasons: 

• The damaged zone where no reduction is performed is small compared to the size of the structure 
in realistic engineering problems; 

• The local fine scale solution is efficiently obtained by Krylov corrections of a good initialisation 
provided by the restriction of the reduced model to the damaged zone. 

Yet, some issues need to be addressed before global/local reduction algorithms can be efficiently 
and safely used to tackle engineering problem: 

• Cost of the construction of the reduced systems. A system reduction is needed to ensure the 
computational efficiency of the method. It is now well known that projection-based model order 
reduction only reduces the size of the algebraic systems to solve, and that the construction of these 
systems is still expensive. Indeed, it involves the evaluation of the internal forces to construct the 
tangent stiffness and the residual of the reduced problem, which requires global updates of the 
constitutive behaviour, and integrations over the whole domain. System approximations permit 
to solve this issue. They basically consist in obtaining a cheap, but reliable approximation of 
the internal forces. In [TS], we used the so-called hyperreduction [15] (also introduced as the 
missing-point approximation by [47 in a different context). We found that adaptive reduced 
basis methods need to be carefully tailored when used within this framework. Indeed, one 
cannot rely on the fact that the residual of the full-scale problem is known at any iteration of 
the nonlinear solver anymore: the number of global evaluation of the internal forces must be 
kept minimal. In addition, we observed that the stability of the particular system approximation 
that we used is highly dependent on the choice of its parameters (reduced integration domain) , 
especially when the damage mechanism localises. We foresee that excluding the most damaged 
zones from the reduced domain can significantly improve the stability of this scheme. Some 
numerical investigations are required to confirm this intuitive speculation. 

• Quality control. Each step in the algorithm is a source of error: (i) choice of the set of degrees 
of freedom for which the solution is searched in a reduced space (ii) choice of the reduced basis 
(iii) choice of the accuracy for the local solution (iv) choice of the system reduction (v) choice of 
the time stepping parameters. These choices should be made within an encompassing framework 
to ensure that the error associated with each of them can be controlled in the same manner. 
Therefore, one needs to find a relevant error criterion and build a control process, capable of 
ensuring that the above choices introduce the same order of error in the solution so that the 
global error be not dominated by any of them in particular. 

In terms of further applications, this method should be validated using more advanced models for 
brittle fracture, in particular nonlocal continuum damage or softening plasticity models, for which 
the distinction between local and global effects is not trivial. The extension to ductile failure and/or 
geometric nonlinearities should also be possible, the global corrections being expected to efficiently 
handle the "smooth" part of the nonlinearities, while the local/global scheme should treat separately 
the localisation of defects. 

We believe that this local/global reduction framework is an important milestone to obtain an 
optimal reduction method applicable to the parallel simulation of failure in solids. Indeed, as alluded 
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to in the introduction of this paper, several authors have proposed to perform a systematic model order 
reduction of the expensive local problems (so-called "problems by substructure") classically involved 
in domain decomposition methods [32 ,33,29. The application of such ideas to failure assessment has 
never been tackled so far, as it first requires an efficient reduction scheme for the local problems, such 
as the one proposed in this paper, in order to guarantee reasonable scalabilities and load balancing. 
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